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Abstract 

We present a semi-numerical algorithm to calculate one-loop virtual corrections to scattering 
amplitudes. The divergences of the loop amplitudes are regulated using dimensional regularization. 
We only treat in detail the case of amplitudes with up to five external legs and massless internal 
lines, although we believe the method to be more generally applicable. Tensor integrals are reduced 
to generalized scalar integrals, which in turn are reduced to a set of known basis integrals using 
recursion relations. The reduction algorithm is modified near exceptional configurations to ensure 
numerical stability. To test the procedure we apply these techniques to one-loop corrections to the 
Higgs to four quark process for which analytic results have recently become available. 
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I. INTRODUCTION 



In order to make reliable theoretical estimates of scattering amplitudes it is important 
to calculate next-to-leading order (NLO) corrections. Such calculations will give a reliable 
estimate of the cross section normalization. In addition, they can also give results for 
kinematic distributions. The latter are important for understanding backgrounds to new 
physics searches in greater detail. 

While most of the two-body final states at hadrori_colliders have been calculated at NLO 

|6| 



and implemented in flexible numerical programs 

mm 

four-body final states are more recent and as yet quite limited 
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ts for three- or 
llQ- One of 



the difficulties is the complexity of the analytic calculations required to evaluate the virtual 
contributions. The computation of leading-order (LO) amplitudes has been automated 



using a variety of methods 



Explicit algorithms for numerical evaluations 



n 



of LO cross sections have been developed |19|, 
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2J| and are of great value 



to experimenters. In contrast, no similar development has occurred for the evaluation of 



NLO amplitudes. Because o 



one-loop integrals 
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recent progres s in the 



recent progres s in the understanding of algorithms to evaluate 

one can envision a similar automation of 



one-loop calculations. This would enable us to extend NLO calculations to more complicated 
final states, such as four-body final states (e.g. the production of two pairs of massive 
quarks), five-body final states (e.g. vector boson plus four jets) and beyond. 

In this paper we report on early steps along this path. Section |n] discusses a general 
approach to one-loop calculations. In section III we outline the numerical implementation 
of an algorithm similar to the one developed in ref. 28| for processes with up to five external 
particles. We use the integration-by-parts method jSJ, |35| to evaluate up to rank-five five- 
point tensor integrals with massless internal lines. In the kinematic regions where a Gram 
determinant is small or vanishes, the basic recursive algorithm discussed in Section III 
breaks down. Section HVl shows how to treat these exceptional regions, and gives an explicit 
algorithm for numerical implementation. In particular, it discusses how to extend the basic 
method, (outlined in ref. 29|), to treat five-point tensor integrals in the exceptional regions. 
The results are shown to be stable close to exceptional phase space points. Section V 
provides a brief outlook. Appendix ^ collects results for the basis integrals. 

The combined algorithms of Sections IIIII and llVI have been implemented in a numerical 
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FIG. 1: The generic A^-point loop graph. 



program. As a first application, we have used the tensor evaluation package to calculate 
the virtual corrections to Higgs plus up to four partons using an effective Lagrangian ap- 



proach 



El- 



II. ANATOMY OF A ONE-LOOP CALCULATION 



The calculation of virtual corrections to a scattering process follows a general pattern. 
Each scattering amplitude Ai consists of a gauge-invariant sum of Nd Feynman diagrams, 

(1) 



M = J^Ai 



i=l 



Each one-loop Feynman diagram A can in turn be decomposed into a sum of contractions 
between two types of tensors. The first tensor is the kinematic tensor K^^...^j^j containing 
couplings and the particle properties such as kinematics, spin, polarization and color. The 
other tensor is a one-loop tensor integral which only depends on the event kine- 

matics, i.e. it depends on the momenta and masses of the scattered particles. The explicit 
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decomposition of each Feynman diagram is given by 

TV 

A{pi,...,Pn;£i,---,£n) = I^'"'''^'{D;qi,. . . ,qN) , (2) 

TV/=0 

where Pi and Si denote the momenta and polarizations of the external particles and the 
tensor integral ' jg given by 

/—W,,.....,.)-/^!^^. + ..^±n, (3) 

and D = 4 — 2e. Note that we are considering the case of massless internal propagators, 
although the method can be extended to include internal masses. The generation of the 
kinematic tensor K, is analogous to the generation of leading-order (LO) scattering am- 
plitudes. That is, the kinematic tensor is composed of LO multi-particle sources. Many 
methods have been developed to deal with the generation of LO amplitudes, either analytic 
or purely algorithmic. Such methods can be adapted to generate the kinematic tensor. 

The more problematic issue is the evaluation of the tensor integral. Using the generic 
expansion of a tensor integral 

TV 



TV 



TV 



+ 2^ aoooois-iiv/fi' 9 Qis Qi^ ■ ■ ■ Qij,! 

+ ■■•, (4) 

we can decompose the tensor integral as known tensor structures and form factor coefficients 

actors, we can evaluate the tensor 



'^hh-ihi- Oiice we have a procedure to calculate the form 
integral numerically using a computer program (see also | 

Using integration-by-parts decomposition techniques each form factor coefficient is iden- 
tified with one generalized scalar integral [4^, as pictured in Fig. ^ 

f I 1 

I{D; {qi, ui}, {qn, i^n}) = I{D\ i^i, V2,...,Vn) = J d"' ■ ■ ■ <f^ ' 

for which a recursive algorithm can be set up in such a manner that all generalized scalar 
integrals can be numerically reduced to the set of analytically known master integrals. Using 



this method we can construct numerical algorithms to calculate the general tensor integral. 
Note that, due to the fact that the master integrals are evaluated using analytic expressions, 
we do not anticipate any significant loss of accuracy for the numerical evaluation of the 
amplitude compared to a purely analytic method. 

One complication is the fact that the evaluation has to be performed within dimensional 
regularization. Because of this we cannot simply calculate a complex valued tensor. Instead 
we calculate the complex valued coefficients q of a Laurent series: 

«iii2-jA/ = ^ + ^ + Co + cie + C(e^) , (6) 



as was advocated in ref. j^. As the tensor integral contains terms of order and 1/e 
we need to specify a consistent scheme for evaluating the kinematic tensor up to potential 
terms of order e^. 

All the elements for a fully automated evaluation of one-loop virtual amplitudes are now 
in place, although they are not yet implemented in a single procedure. Indeed, the method 
we used relies on the generation of the graphs using Qgraf and the simplification of 
terms proportional to the square of the loop-momentum using the algebraic manipulation 
program Form This latter step is a necessary prerequisite to evaluate the kinematic 
tensor. Were terms proportional to g'^'^ to remain in the kinematic tensor, they would lead 
to ambiguities in performing the contraction with the tensor integral. ^ 

For amplitudes involving more than one quark line one can encounter terms of the form 
u{pi) .../... u{p2) u{j>z) •••/••• uipi), (where / = •jf^l^), which can not be simplified. For 
the simple case H qqq'q' (and H — > qqqq) we performed the square analytically, thus 
removing all spinor lines. After elimination of the spinors, the dangerous terms quadratic 
in / can always be removed. 



III. NUMERICAL IMPLEMENTATION OF THE ALGORITHM 

In this section we will outline the numerical implementation of the evaluation of tensor 
integrals fi^'a ■ma/j^/}- q2, . . . , qj^) in terms of a Laurent series. The first step is to reduce 
the form factor coefficients a,j^j2---jM generalized scalar integrals /(-D; i^i, z/2, . . . , vn)- 

^ By eliminating all metric tensors from the kinematic matrix we ensure that the kinematic tensor is 
composed entirely of external vectors which are purely four dimensional. 
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Next, the generalized scalar integrals have to be reduced to the set of master integrals. 
Finally, the master integrals have to be evaluated. 



A. The tensor integral decomposition 



The tensor integral decomposition is performed using the method of Davydychev |4ll |. 
which is derived using the integration-by-parts method. We introduce the notation for the 
generalized scalar integrals, eq. (jSl), 

I{D;{iyk}ti) = HD;i^i,i^2,...,i^N). (7) 

The general form for the tensor integrals with external legs and M free indices can be 
expressed in terms of these scalar integrals 

i,,...,AD;{Qu'^i}---,{QN,'yN}) = (-0 {[9]'[qir---[qNr^},....,M (8) 

A.ftl,K2v:'^iV>0 

X . . . {i^n)kn + 2(M - A); vi + ki, . . . .vn + i^n) , 

where 

M. = ^. (9) 

is the Pochhammer symbol. The structure • • • [(lN]'^'^}fii...fii^i is a tensor which is 

completely symmetric in the indices fii, . . . , fiM and constructed from A metric tensors g, ki 
momenta pi, K2 momenta p2, etc. Once we have a method to evaluate the generalized scalar 
integrals, eq. (jHl) can be used to construct the tensor integral. 

Thus, for example, the rank-1 tensor integral with N external legs can be decomposed as 

TV 
ii=l 

= + 2; 2,1, !,...,!) +/(/} + 2; 1,2, !,...,!) 

+ ■■■ + /p + 2;l,l,l,...,2) (10) 
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Similarly, a rank-2 tensor integral is given by 

2! N 

= -^J(D + 2;1,1,1,...,1) (7^^^^ 
+2 /(Z} + 4;3,1,1,...,1) gfigf 
+ + 4; 2,2,1,...,!) (g^gf +grg2^^) 
+ •■•, (11) 

where the first sum runs over the permutations of the indices P{yUi,/i2}. 
B. Recursion scheme for scalar integrals 

We will outline and improve the recursive scheme proposed in ref. [2^ , presenting explicit 
formulas for the case of five or less external particles. Using the basic equation of the 
integration-by-parts method, an identity which is valid for arbitrary values of the parameters 
Hi can be derived. 



0. (12) 



Differentiating we obtain the base identity 



N / N \ 

Y.[Y.S,.yAu,I{D-{u, + 5k,}ti) 
j=i \i=i ) 

N / N \ / N 

- J2 ydiD - 2; {u, - 4.}f=i) -[D-i-J2''A[Y.y^] ^(^; Mk=i) ^ (13) 
1=1 \ j=i J \i=i J 

where S" is a kinematic matrix which, for massless internal particles, takes the form 

Sij = {qi - qjf ■ (14) 

The dimensional shift identity 

N 

I{D - 2; {u,}ti) = -J2 ^^^(^-^ + ^^'^}^=i) ' (15) 



i=l 



2 3 4 5 6 

D/2 



FIG. 2: Reduction scheme for triangle diagrams with three off-shell legs. 

was used to obtain the base identity. 

The integrals are characterized by the value of n, which is the naive degree of ultraviolet 
divergence 

(16) 



D 

n = I a 

2 



e=0 



It is also helpful to introduce some further notation 

N N N N 

^-E^' ^^-E^'' ^-E^^ = E^'- (17) 

1=1 j=l i=l *ii=l 

The basic recursion relation is obtained by solving SjiUi = 6ij in eq. (fT^ and shifting ui 
by one unit. The result is, 

{u, - 1)I{D; 

N 
i=l 

Figs. 121 El and m illustrate the action of the recursion relations in the plane defined hy D/2 
andcr. The location of the basis integrals is indicated by a diamond. The recursion relation 
Eq. ()18p either reduces both the dimension D and the value of a by 2, or it keeps D fixed 
and reduces cr by 1. Its action is indicated by the red (dotted) lines in Figs. El and El 

By performing the shift {uk\k=i ~^ {^k + Sik}^=i, summing eq. (fT8|) over the index /, and 
using eq. (fTK|l we derive the standard recursion relation [4^, 

N 

(D-l-a) BI{D- {z/Jti) = - 2; {z/.}f=,) -Y.hJ{D - 2; {v, - 5,,}^^. (19) 

1=1 



Eq. ()19p reduces D by 2 and cr by or 1. Its action is illustrated by the blue (solid) lines 
shown in Figs. |21 El and El 

For the case of pentagon diagrams with all = 1 and D = 4 — 2e, we will use the 
following recursion relation to arrive at the basis pentagon m D = Q — 2e dimensions, plus 
boxes 

TV 

I{D- {u,}t,) = {D + l-a) BI{D + 2- {u,}t,) + ^ ^^^(^^ " ^^^^k=i)- (20) 

i=l 

Its action is illustrated by the magenta (solid grey) line in Fig. jH 

By combining eqs. ()18j) and (fTn|l we can eliminate the term with a prefactor depending 
on the dimension D, to yield a fourth recursion relation, 

{u, - 1)I{D; 

= -|/(^ - 2; {u, - 6t,}t^) + E - SrA I{D - 2; {u, - 6,, - 6t,}ti). 

i=i ^ ^ 

(21) 

This identity is particularly efficient since it reduces D by 2 and cr by either 1 or 2 units. 
Its action is indicated in Figs. El El and jH by the green (dashed) lines. 

Note that for the case = 3 in Fig. |21we assume that det(5') ^ 0, which is only the case 
if all three external momenta are off-shell. Triangles with only one off-shell leg are basis 
integrals. The case of triangle integrals with two off-shell legs will be treated in the following 
section. 



C. Reduction of two mass triangles 



For the specific case of two mass triangles (pi,p| 7^ 0,p| = 0), the matrix 



S 



( vl 
v\ 

\v\^ y 



(22) 
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D/2 



2 3 4 5 6 

D/2 



FIG. 3: Reduction scheme for box di- 
agrams. 



FIG. 4: Reduction scheme for pen- 
tagon diagrams. 



is singular. Nevertheless useful recursion relations can still be derived from eq. |T!?|) . By 
choosing Ui which satisfy the relation, Yli SjiVi = CLj 







( " ] 







, y = 


a 

pI 


(23) 


vv 




\1^J 





and choosing the parameter a such that Ui = 0, we derive the recursion relation valid 
for Ui > 1 and pf ^ p\ 

I{D- z/i, 1/2, i^s) = 2 ^ 2/ \a i^i-l. V2-\ v^)-I{D-2- 1/2, 1^3-1)1 ■ (24) 

This relation lowers -D/2 by one unit and o by two units and is shown by a green (dashed) 
line in Fig. El 



By making the choice for a, 



/ n \ 



(25) 



we derive a further relation 

1 1 



i/i, z/2, z/a) 



p\ liD - 2; z/i, 1/2 - 1, z/s) - P2 -^(-D - 2; z/i, 1/2, 1/3 - 1) 



(26) 
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2 3 4 5 6 

D/2 



FIG. 5: Reduction scheme for triangles with two ofF-shell legs. 

valid for p\ p\. This relation lowers D/2hy one and a by one and is shown by a blue 
(solid) line in Fig. El 

Note that there is a potentially dangerous case when vi = 1 and = a + 1 — 2e, however 
if we restrict our attention to = or n = — 1, this possibility can only occur for D = 4 — 2e 
and a = 3 which is a basis integral. 

Finally we note that in the example considered later, Higgs decay to four partons, p\ = p\ 
can never occur for two-mass triangles. The degenerate case pi = p\ can however also be 
treated easily if needed. For instance, by solving SijUj = aj with J2j Vj = 









f \ 


a = 





, y = 


a 




VJ 




[-a) 



one obtains 

I{D; z/i, z/2, z/3) = -I{D; + 1, 1.3 - 1) . (28) 
This relation can be used to reduce the triangle to generalized self-energies. 

D. The standard algorithm 

We summarize here the detailed steps of the algorithm as we implemented it. In many 
cases there is freedom in the choice of the particular recursion relation used. 
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The set of master integrals, which are the end-points of the recursion relations, consists 

• all tadpoles I{D\vi) for arbitrary D and vi ; 

• all two-point functions /(-D; i^i, 1^2) for arbitrary D and {i^k} 

• all one- mass triangles i^i, z/2, z/3) for arbitrary D and {uk} ] 

• two-mass triangles I{D = 4 — 2e; 1, 1, 1) ; 

• three-mass triangles I{D = 4 — 2e; 1, 1, 1) ; 

• all boxes I{D = A- 2e; 1, 1, 1, 1) ; 

• all pentagons I{D = 6 — 2e; 1, 1, 1, 1, 1) . 

For any given generalized scalar integral the algorithm proceeds as follows: 

• check if the integral has already been computed, if so read in the value of the integral, 
which has been previously stored; 

• check if the integral is a basis integral, if so evaluate its complex Laurent expansion 
using the analytic results (see Appendix and store the result; 

• if iV = 3 and det(S') = apply eq. (pUj) (reduction of two- mass triangles); 

• if a = AT = 5 and D = 4 - 2e apply eq. ((201); 

• ifcr-iV + 2- e = f apply eq. (HHD; 

• if a-iV-F3-e = f and o->iV apply eq. dH)); 

• ifa-A^ + 3- e = fanda- = A^ apply eq. (UHl); 

• ifa-iV + 4- e<f apply eq. 

This procedure is repeated until all generalized scalar integrals have been evaluated. 
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IV. THE MODIFIED ALGORITHM FOR EXCEPTIONAL KINEMATIC RE- 
GIONS 



The numerical implementation of the recursion algorithm as outlined in the previous 
section assumes that the kinematic matrix S is invertible for = 4, 5 and the quantity B is 
not zero. If either B or det(S') is close to zero the numerical algorithm of section III becomes 
unstable. 

For these exceptional kinematic regions of phase space we need to develop numerically 
sound procedures. Such exceptional regions are usually associated with thresholds (i.e. a 
massive particle produced at rest) or planar events (i.e. a linear dependence between the 
particle vectors). 

The exceptional kinematic regions are characterized by a vanishing Gram determinant. 
The Gram matrix is defined as (assuming = 0) 

G,,=2q,-qf, it,j = l,2,...,N-l) , (29) 

and is closely related to the kinematic matrix 

Sij = Wi ~ IjY — ~'~ 2^-'"' ~ ~ ^ ^ ' ^^^^ 

The determinants of S and G are related through the B quantity of eq. ()17|1 

det(G) = (-1)^-^5 det(5) . (31) 

where det(G) denotes the determinant of the (A^ — 1) x (A^ — 1) submatrix. We introduce 
some arbitrary small cutoffs Go, 5*0. We then distinguish 4 cases: 

• Both det(5) > 5*0 and det(G) > Gq: 

In this case we apply the standard recursion relations of section III. 

• The det(^) < and det(G) > Go region: 

This region is in fact not a kinematic exceptional region. We will rewrite the recursion 
relation in terms of the invertible G matrix. This gives an alternative set of recursion 
relations which do not depend on the kinematic matrix S. 

• The det(S') > 5*0 and det(G) < Go region: 

In this case the inverse of S is still defined, while the B parameter becomes small. 
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we al- 



This renders the recursion relations eqs. p9|) and (j2U) unusable. In ref. 
ready outlined a method to deal with this situation by rewriting the affected recursion 
relations as an expansion in the small parameter B. 

• The region where both det(S') < 5*0 and det(G') < Gq: 

In this region all recursion relations of section III have to be rewritten as expansions 
in the smallest eigenvalue of the kinematic matrix S. 

A. The region det(S') < 5o and det(G) > Gq 

The numerical procedure of section III becomes unstable because the inverse of the kine- 
matic matrix S is ill defined. We encountered exactly this condition for the double off-shell 
triangles. As we saw in that case two types of equations can be derived which we need to 
generalize to boxes and pentagons. 

The first set of equations is obtained by looking for the eigenvectors associated with 
the smallest eigenvalues of S*. If the determinant is exactly zero one can find a non trivial 
eigenvector y in the null-space of Sij and use the base identity eq. (fT^ to derive 



N / ^ \ 

= - 5^ yJiD - 2; {u, - 6.,}ti) - (/^ - 1 - a) ■ (32) 

i=l \i=l J 



If ^^j^ 7^ one obtains the relation 

i=i z^i=i y 



1 ^ 

m Mti) = E ^^^^(^ - 2; {^k - 5..}f=i) , (33) 



which reduces both D and a (and possibly A^), while keeping n fixed. 



If J2iLi Z/« = one obtains 



i{D- Mt,) = -Yl -m - + 5.af=i) . (34) 

This relation is not very efficient since it reduces neither D nor a. However repeated use of 
eq. with I fixed, is guaranteed to reduce A^. 

We now discuss how to extend eqs. (jH^ and to cases where the determinant is small, 
but non zero. In this case one computes the eigenvector y of the matrix Sij which has the 
smallest eigenvalue. The terms 

N / N \ 

j=l \i=l J 
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in eq. p3p give only a small contribution to the equation. This means that, in the case 
Sill 2/j 7^ 0, one can rewrite eq. (fT^ as 

1 ^ 

^ - 1 - ^ Tit E.=i 

+ (36) 



This equation allows one to reduce the integral on the l.h.s. to a sum of "simpler integrals" 
(the first term in the r.h.s.), i. e. integrals with lower lower a and possibly lower A^, 
plus a sum of more "complicated integrals" (the second term in the r.h.s.), which have the 
same D and higher cr, but which are suppressed by the small coefficients SjiUi. Notice 
that since these integrals have exactly the same kinematic matrix as the starting integral, 
the computation of these integrals will make use of the same relation eq. (jHUI) . with only 
modified D and exponents {i'k}k=ii giving again simpler integrals and "more complicated 
integrals" which are suppressed by the square of the small coefficients SjiUi. At some 
point the small correction terms will be less than the desired accuracy and the reduction 
can be terminated. 

Similarly, in the case Yl!i=i Vi close to zero, one obtains 



- (D-a) ^^^I{D + 2; {u, + 5i,}ti) , (37) 

where now both the second and the third term in the r.h.s. are small corrections to the 
result. For both integrals the kinematic matrix is again the same as the starting one, so 
that a repeated application of this same identity allows one to achieve the desired accuracy 
in the answer. 

Notice, that if D = 1 + o" — 2e in eq. then the prefactor is of O (1/e). Therefore to 
have a correct result for I{D; {z/^}^^) up to O (1) one needs the e-expansion of the integrals 
in the r.h.s. up to O (e). 

To avoid O (e) expansions in the case D = 1 + a — 2e we derive an alternative set of 
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recursion equations. Using eq. ()3U|) we write 

N 1^1 ^^"^ 

j=i j=i j=i j=i 

where we parametrized the scalar integral such that = 0, i.e. GiN = 0. This means we 
are free to determine the values of {i = 1,2, . . . , N — 1) by solving Xljl^^ ^ijVj — 
any non-zero vector r of our choosing and to define Un = — Yl!j=i Vj- have 

N 

Y,Si,yj = Cir)x{l,l,...,l)-r, (39) 

where C(r) = | J^fJi GjjVj = Ef=~/ QjVj = J2Z'=i ^i^^^^j 
The new base equation now becomes (with g^v = 0) 

N-1 N-1 
j=l i,j=l 

-C{r)I{D - 2; (40) 

To derive usable recursion relations we chose rp'' = 5ik {k = 1, . . . , N — 1) and find the 
set of (A^ — 1) recursion relations 

Af-l 

(z/fc - l)m {ui}l,) = Y G;,' {I{D - 2; {ui - 5u - 5i,}l,) - I{D - 2; {ui - - 5iu}l,)) 

i 

-CkHD - 2; {ui - 6ik}l,) (41) 

withC, = j:l'i'QfG-,\ 

By summing over k we derive the recursion relation for the case that uj^i ^ 1 

N-1 

(z/AT - 1)I{D- {ui}^,) = - Yg^ {I{D - 2; {ui - Su - 5,^}£i) - I{D - 2; {ui - 26in}1,)) 

i 

+ {C-l)I{D-2-{vi-5ir,}l,) (42) 

with C = J2k=i Ck and gi = J2k=i G^k- 

Note that these equations can be applied only if at least one of the z/j > 1. The cases 
where all Vi = 1 and D = 1 + cr — 2e correspond always to basis integrals, so that all cases 
are covered by the recursion relations presented here. ^ 



^ In our current algorithm wc chose to apply cq. (|36() for triangles even \i D = 1 + a — 2e. In this case we 
need the expansion of two-point functions up to O (e). Suitable results are reported in Appendix IXI 
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B. The region det{S) > Sq and det(G) < Go 

In the case where B vanishes but det{Sij) ^ problems arise only if one is to use relations 
eqs. (jl9|) or (j2H) . while eqs. (jl8j) or (j2(J|) can be used in a straightforward way. 
For small B one can simply use eq. (|2(jp 

N 

I{D- {vk}Li) = E ^^^(^^ {^'^ - ^^'^>^=i) + {D + l-<y)BI{D + 2- {z/Jt,) , (43) 

i=l 

where the first term gives rise to simpler integrals (lower a and possibly lower A^), while the 
second term gives a small correction to the result and is computed by applying this same 
recursion relation iteratively until the desired accuracy is reached. 



C. The region det(5) < 5*0 and det(G) < Gq 

First we consider phase space points where the Gram determinant exactly vanishes. In 
this case we have a non-trivial solution of the equation 

7V-1 

G.,z, = . (44) 

Setting zn = — Yl!i=i ^* leads to 

N ^ N-l Af-1 

s^,z^ = c X (1, 1, . . . , 1) , c^-Y. ^^^''^ = E ^^'^'^ • (^^5) 

j=l k=l k=l 

If C > Co we can rescale zf. Vi = Zi/C. From eqs. (fT^ and (jlSI)) we obtain 

N 

■t=i 

If det(Gjj) is close to zero, the equation SijVj = 1 will be approximate, leaving small 
correction terms A, 

5',,r, = (l, !,...,!) + A,. (47) 

By explicit construction it is possible to keep exact the condition XljLi = 0. In this case 
the expanded recursion relation becomes 

N 

1=1 

N 

+ Y ^^'^^HD + 2; {uk + 6,k}Li) ■ (48) 



1=1 
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If C < Co we cannot rescale the vector Zi. Instead we write an expansion 

- E ''"^ (^^- + ^^■')^(^ + 2; {^^- + + ^'4f=i) • (49) 

■ 1 ' 

We finally note that the expanded recursion relations change the value of n defined in 
eq. (fTH|l . Since more steps might be needed to reach the desired accuracy these relations can 
lead outside the normal IR and UV boundaries, giving rise to deep-IR or deep-UV integrals. 



D. Summary of the modified algorithm 

The set of standard and expanded recursion relations is over-complete. Therefore there 
is some freedom in most steps of the algorithm. We summarize here the detailed steps of 
the modified algorithm, as we implemented it. While the procedure works, it is yet to be 
optimized. 

The set of master integrals is the same as the one listed in section UlI Dl Notice that, as 
explained above, two-point functions need to be expanded up to O (e) and basis integrals 
out of the standard IR/UV boundaries will appear. 

We rescale the hard event by a typical hard scale of the process. In order to measure 
the closeness to the exceptional momentum configuration we introduce a small parameter 
C ^ 1 (which we fixed at C = 10~^ in the numerical examples). 

Given any scalar integral, the modified algorithm proceeds as follows: 

• check if the integral has already been computed; if so, read in the value of the integral, 
which has been previously stored; 

• check if the integral is a basis integral; if so, evaluate its complex Laurent expansion 
and store the result; 

• if all eigenvalues of Sij are larger than ( the standard algorithm of section IIII Dl can 
be applied, with the following modifications: 

. ifa = A^ = 5andD = 4-2e, or(T-A^ + 2- e>f apply eq. 
■ ifa-iV + 2- e= f apply eq. (UHl); 
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■ ifa-iV + 3- e= fandB>C apply eq. (EH); 

■ if (ct - + 4 - e < f or (a - iV + 3 - e = f and a = A^)) and 5 > C apply 
eq. (HI; 

■ ifa-iV + 3- e= f and a = N and B > C apply eq. (O; 

■ ifa-A^ + 4- e<fand5>C apply eq. (HHl); 

■ ifa-Ar + 3- e<fandB<C apply eq. PH) . 

if some eigenvalue of Sij is smaller than C, but all eigenvalues of Gij are larger, the 
event is not really exceptional, however standard relations of sec. IIIIDI can not be 
applied since the kinematic matrix has no inverse, accordingly 

■ if D 7^ a + 1 - 2e or iV = 3 apply eq. (jSEI); 

. if D = o- + 1 - 2e, iV 7^ 3 and some z/fc 7^ 1 (A; = 1, . . . , iV - 1) apply eq. (gl}; 

■ if = a + 1 - 2e, iV ^ 3 and I/7V ^ 1 apply eq. (j42|) . 



• if both Sij and Gij have eigenvalues smaller than C,: the event is considered exceptional. 
In this case 

■ if C > apply eq. (jlHj) (where C is defined in eq. (jl^ ): 

■ if C < apply eq. P^. 

The number of terms included in the expanded relations is limited by the requirement that 
^mtcr-i ^ Here niter denotes the number of iterations and ki are the small coefficients in 
eqs. (f^UI) . (jHUj) . (PHj) and (jlHI), which parameterize the departure from the exceptional phase 
space points. Concretely, the ki are the prefactors of the integrals I in the second term in 
the r.h.s of those equations. 

E. Some explicit results for H — > qqq'q' 

An ideal process to illustrate expansions close to exceptional momentum configurations 
is the decay 

H{-P5) q{pi) q{P2) q'ip-i) q'{Pi) ■ (50) 



19 



This is a non-trivial example, with N = 5 and one massive external momentum, where 



the analytical result has recently become aval 



process H qqqq which is also given in ref. 



able 



• 3)- 



(Alternatively one could study the 
A close inspection of this result reveals 



immediately that it is free of any spurious singularities and therefore its evaluation is stable 
for any phase space point. ^ 

We present here a comparison between the numerical and the analytical result for phase 
space points approaching exceptional momentum configurations. 

For the Higgs decay process we are considering, the determinant of the kinematic matrix 
defined in eq. (IT^ is given by 

det(S'ij) = 2Si2S23S34(Sl23S234 " pIs2z) , (5l) 

where Sijk = Sij + sjk + Sik and the momentum assignments are given in eq. (jMlj) . It is 
straightforward to see that if is a linear combination of pi and p2 + Ps, the determinant 
of the kinematic matrix vanishes. 

Pi = api + [3p2 + 7P3 with a = — — and 7 = /3 =^ detfS'j,) = , (52) 

psi2 + 7S13 

where the condition on a ensures that p4 is massless. It is then possible to approach the 
limit det(5'jj) = slowly, by choosing 7 = /? + 5, where 5 is a small variation such that 
5 <tia,(3, 7. 

As an illustration we present results for the following exceptional momentum configura- 
tions, {E,p^,py,p^), 

= (-0.50000000000,0.00000000000,0.00000000000,-0.50000000000) , 
P2 = (-0.50000000000, 0.00000000000, 0.00000000000, +0.50000000000) , 
P3 = (+0.10000000000, 0.05000000000, 0.00000000000, +0.08660254037) , (53) 

Pi = api + pp2 + 'JP3 a = -p'jS23/iPsi2 + 'jSi3) /?=l/3 7 = /3 + (5, 

P5 = -Pl -P2-P3-P4, 

where (3 is an arbitrary number O (1) and we vary S in the range —50 < ln5 < 0. 
In Fig. El we plot the absolute value of the relative accuracy 

p=^Mf^, (54) 

Av,A 



Remember we are interested in separate jet-production, i. e. we always assume \sij \ > Smin > 0. 
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Ln{det{Sij)) 



Ln(det(Sij)) 



30 -25 -20 -15 

Ln(det{Sii)) 



FIG. 6: Relative accuracy \p\ for the pole (left), the 1/e pole (center) and the constant part 
(right) of the one-loop amplitude squared for H — > qqq'q' as a function of the determinant of the 
kinematics matrix Sij. 

where, as defined in 0], Av,n and Av,a denote the numerical and analytical result for the 
full one-loop amplitude squared for the process H qqq'q' . 

The matrix Sij denotes here the kinematic matrix for the standard ordering of momenta, 
as defined in eq. (fT^ . In the evaluation of a full one-loop amplitude one encounters {N — 1)! 
such matrices, which are obtained by permuting (A^ — 1) momenta keeping one arbitrary 
momentum fixed 

S^J = - qa{j)f . (55) 

where a{i) denotes the permutation of index i. Additionally, external momenta are con- 
tracted together, giving reduced matrices. Whenever any such matrix is close to singular, 
the modified recursion relations are used as explained in sec. IIVDI Therefore plotting the 
accuracy as a function of just one determinant is only indicative. 

From Fig. El we see that the numerical answer is well-behaved when approaching excep- 
tional phase space points, det(S'jj) 0. Since the closeness to the exceptional phase space 
point is the expansion parameter in the modified recursion relations, an extremely good 
accuracy can be reached very close to singular points by including only a few terms in the 
expansion. In the intermediate region, where the expansion parameter is larger, one has to 
include more terms in the expansion. This however means computing a larger number of 
integrals with higher D and a. In the current example, eq. ()50|) a relative accuracy of 10~^ 
was achieved, limited by the computer memory. A greater accuracy could be obtained by 
allocating more computer memory. 

We now study the case i? — > 0. If the four- vector is a linear combination of pi,p2 and 
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Ln(IBI) 
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Ln(IBI) 




-15 -10 
Ln(IBI) 



FIG. 7: Relative accuracy \p\ for the pole (left), the 1/e pole (center) and the constant part 
(right) of the one- loop amplitude squared for H — > qqq'q' as a function of ln(|i3|) for the box 
integral obtained by setting 1/3 = 0. 

P3 then for any choice of a, (3 and 7 which keeps pi = 0, we find B = 0. However for the 
process in eq. ()50|) the only pentagon integrals appearing in the numerical answer are scalar 
integrals in D = 4 — 2e with cr = 5. These integrals can be reduced using the standard 
eq. (pni) which holds regardless of whether i? = or not. 

To study the B expansion we therefore choose specific phase space points close to 7 = —1. 
In this case the coefficient B of the reduced 4x4 kinematic matrix, which is obtained by 
setting z/3 = 0, vanishes and an expansion in B is needed. We choose the same momenta as 
in eq. (j^Hj) . but fix 7 = —1 + 5 and vary 6 in the range —8 < In 5 < 0. 

Fig. [7| shows then the absolute value of the relative accuracy for the numerical one-loop 
amplitude squared for H qqq'q' , plotted as a function of ln(i?) of the reduced matrix 
obtained by contracting the momenta p^. and pi (z/3 = 0). Again, numerical results are 
stable close to exceptional momentum configurations and a relative accuracy of 10~^ can be 
achieved. 



V. OUTLOOK 



We presented here a semi-numerical method to calculate one-loop corrections to scattering 
amplitudes with up to five external legs and massless internal lines. A variety of approaches 
to deal with the exceptional kinematic regions have been developed. Further developments 
are possible, guided by practical applications. We demonstrated that numerical methods 
can be extended to kinematic regions such as threshold and planar events. 

With these results we can integrate the virtual matrix element of js^ over all of phase 
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space to obtain the fully differential Higgs boson plus two jet cross section at next-to- leading 
order. 

It is now straightforward to extend the calculation to other, yet uncalculated, cross 
sections. Examples of processes of current experimental interest are di-boson plus one jet 
{Vi, V2, j), tri-boson production (Vi, V2, V3) and vector boson plus heavy quark pairs (VQQ). 
The extension to the case with internal masses can be implemented without much difficulty. 

Furthermore, the recursive techniques are not restricted to 5 external particles. An 
extension of the method to 6 or more external particles is only limited by computer resources. 
To what extent this will be feasible in practical applications still needs to be investigated. 

Finally, because the entire calculation is done numerically we have great flexibility in 
recalculating the cross sections in any "scheme" necessary. For example, this flexibility 
should be helpful in combining NLO calculations with shower Monte Carlo programs. 

APPENDIX A: INTEGRALS 

We report here expressions for all of the basis integrals needed as end points for the 
recursion relations used in our approach. All of these integrals can be found in the literature, 
but we find it useful to collect the results here. 

We first introduce some notations. We will denote In{D; ui . . . v^) the master integral 
/(-D; vi . . . un) with N external legs. The quantity cr which is equal to unity in 4 dimensions 
is common to all integrals and is given by 




(Al) 



We will frequently need the following expansions in e 



r{n - e) 
r{-n-e) 

r(i-e) 




i^l - eH.in + 1) - Y {H^n + 1) + H2{n + 1)} 



(A2) 



where the generalized harmonic sum is defined by 




n-l 



(A3) 



m=l 
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FIG. 8: The generic self energy diagrams. 
The dilogarithm is defined as usual as 

r df 

U^{x)^- / ^ In(l-t). 
Jo ^ 

1. Tadpole diagrams: Ii{D = 2{cj + n) — 2e; vi) 

Within dimensional regularization tadpole integrals are zero 

d^l 1 



(A4) 



r d^i 1 

h{D;.,) = J -1571^ = 0. (A5) 



2. Self energy diagrams: hiD = 2{a + n) — 2e; J^i, 1^2 I {p?}) 

All generalized scalar two-point integrals, shown in fig. |Hl are master integrals. 

2.. r d^i 1 



hiD = 2(c7 + n-e);z/i,i/2| {pf}) 



%1T 



o-,n-e^i.^i + ^ - e)r(z/2 + n- e)V{-n + e) 



= i-lYi-pir-' ' ' \ . ' \ / \ (A6) 

^ ^ ^ r(z/i)r(z/2)r(a + 2n - 2e) ^ ^ 

Depending on the value of n, defined in eq. ()16|) . we have different expansions. The cases 
n > and n < — 1 are needed to deal with the exceptional momentum configurations. 
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The UV case n > 0: 

l2{D = 2{a + n-e);u,,U2\{pl}) 
1 



[-ly+^cr r(z/i + n)T{u2 + n) {-pj) 



r(z/i)r(z/2)r(n + l)r(a + 2n) 
X < - + 2ifi(a + 2n) + Hi{n + 1) - Hi{ui + n) - Hi{u2 + n) - \n{-pl) 



+ 



2Hi{a + 2n) + + 1) - + i^i) - + z/2) - ln(-p2) 



4i72(c^ + 2n) + if2(?^ + 1) - ^2(n + ui) - H2{n + U2) 



(AT) 



The IR cases n < 



1. 1^1 > — n, z/2 > — n: 



2\" 



/2p = 2(a + n-e);z/i,z/2|{p?}) 



;-l)'^cr r(-n)r(z.i + n)r(z.2 + n) {-pD 



X |l + e 



r(z/i)r(z/2)r(a + 2n) 
2ifi(a + 2n) + ifi(-n) - H^{i^i + n) - Hi{u2 + n) - \n{-pl) 



(A8) 



2. i/i > — n, V2 < ■~''^! z^i + 2/2 > ~2n : 



l2{D = 2{a + n-e)-Vi,V2\{p\]) 



2\n 



(-l)''i+"+icr r(z/i + n)T{-n) {-pD 

r(z/i)r(z/2)r(i - n - z/2)r(a + 2n) 



X <^ - + 2i7i(a + 2n) + - i7i(z/i + n) - i/i(l - z/2 - n) - \n{-p\) 



+ - 



2i7i(a + 2n) + H^{-n) - Hi{ui + n) - Hi{l - U2 - n) - \n{-pl) 



+ AH2{a + 2n) - H2{-n) - H2{ui + n) + H2{1 - z/2 - n) 



(A9) 
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3. z/i > — n, z/2 < — n, z/i + z/2 < — 2n : 



2\n 



l2{D = 2{(T + n-e);ui,U2\{pl}) 



2(_l)'^2+«cr r(i^i + n)T{-n)T{l - a- 2n){-pi) 



+0 {e') 



r(z/i)r(z/2)r(i - n - z/2) 

2Hi{l -(T-2n) + Hi{-n) - Hi{ui + n) - Hi{l - z/2 - n) - \ii{-p^] 



(AlO) 



4. z/i < —n, V2 < —n : 



, r 2(-l)'"+^Cr r(l - a - 2n)r(-n) , 

X < - + 2i7i(l - a - 2n) + i/i(-r2) - i/i(l - n - Ui) - Hi{l - U2 - n) - \n{-pl) 



+- (2iJi(l -2n-a) + Hi{~n) - Hi{l -v^-n)- Hi{l - U2 - n) ~ H-pl)f 

- 4i/2(l -2n-cj)- H2{-n) + H2{1 - n - Ui) + H2{1 - n - z/2) 

+0 (e^) . (All) 



3. Triangle diagrams: I3{D = 2(cj + n - e); z/i, z/2, 1/3 | {Pi,P2'P3}) 

In the case of the 3-point general scalar integrals, shown in fig. El we distinguish three 
cases depending on the number of off-shell legs. 



a. One off-shell leg 

The first case is the three-point generalized scalar integral with one off-shell leg. This is 
again a master integral. 
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FIG. 9: The generic triangle diagram. 
The general expression for the triangle with one off-shell leg {pf = = 0) is given by 



Is{D = 2{a + n)- 2e; i^i, u^. i^sljO, 0,^^}) 



2x„-er(i^2 + n-e)r(i/3 + n-e) r(-n + e) 



cr 



r(z/2) r(z/3) r(a + 2n-2e) 



{a + 2n-l)\ \ r(l + e) 



r(l - 2e)r((T + 2n) 
r((7 + 2n-2e) 



(A12) 



Depending on the value of n, the integral is IR {n = —1), UV (n = 0) or deep UV (n > 1) 
divergent. We will look at these three specific cases: 

• The IR case: n = —1 

The generic integral now becomes 



2\-e 



li^'iD = 2{a-l- e); i/i, i^^) = -cr 



X 



1 



(-1)'^ i-pi) 

r(l - 2e)r(a- 2) 



Vr(z/2)r(l -e)J \T{u^)r{l - e); V r(a - 2 - 2e) 
Depending on the values of z/2 and we can distinguish three cases 



(A13) 



1. 1/2 = 1^3 = 1: 0" = '^i + 2 
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This is the most singular case (C(l/e^)) 



/^^(Z} = 2(1.1 + l-e);z/i,l,l) = -cr 



^,+1{2HM)-H-pI)) 



+ - \n\-pl) - 2HM) H-pI) + 2 {Hl{v,) + H^{v,)) 



+ 0{e) 



(A14) 



2. (z/2 = z/ > 1; 1^3 = 1) or (1^3 = z/ > 1; z/2 = 1): a = z^i + + 1 
This is the less singular case (C(l/e)) 

(-1)^1+^ 1 1 
(^,^ + ^,_2)!(zy-l)^ 



X 



- In(-p2) ^H^{u-l) + 2Hi{iyi + u - 1) 



0{t 



(A15) 



3. 1^2 > 1; 1^3 > 1: O" = Z/i + z/2 + Z.3 
This is the finite case ((9(1)) 



jf^(D = 2(a-l-e);z/i,z/2,z/3) 



cr(-l) 



PK^^I + z/2 + Z/3 - 3)!(z/2 - l)(z/3 - 1) 

+0{e) . (A16) 



The UV case: n > 

In this case we simply get from eq. ()A12p 

jL^(D = 2(a-n-e);z/i,z/2,z/3) = Cr{-lY^^'''^''\plr 



n!(z/i + z/2 + z/3 + 2n-l)! 



X 



In(-p^) ^Hi{n^-l) + 2Hi{vi + z/2 + z/3 + 2n) 



(A17) 



b. Two off-shell legs 



h{D = 4 - 2e; 1, 1, l\{plpl 0}) = ( (-p^) 



2\~e 



2\-~e 



-Pi) 
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P3 



1+qi 




Pi 



i+q4 



FIG. 10: The generic box diagram. 



P4 



c. Three off-shell legs 



The basis integral is the finite four dimensional triangle |4fil . |4 



h{D = A-2e;l,l,l\{pipipi}) 



Li2(-px) + Li2(-p?/) +\n{px)\n{py) 

(A19) 



where x = p\/p% y = pI/pI and 

Kx,y) = {l-x-yf -Axy, p{x,y) 



(A20) 



l-x-y+ ^/X{x,y) 

Eq. ()A19|) is valid in the region where all < and X{x, y) > 0. Details about the analytical 
continuation to other regions can be found in reference j48|. 



4. Box Integrals: 1^(0 = 4 - 2e; 1, 1, 1, 1) 



All 4 — 2e dimensional box integrals with a = 4 belong in the basis set. The basic diagram 



is shown in Fig. ^1 The integrals are defined in terms of Si2 = {pi + P2] 



q^Y and 



S23 = {p2 + PzY = {qs — and the off-shellness of the four legs, pf. 

We give here the explicit expression for the 4-dimensional boxes for massless internal lines 
and any number of external masses. 
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a. Four on-shell massless legs 



h{D = 4 - 2e; 1, 1, 1, 1, {si2, S23, 0, 0, 0, 0}) 



cr 



S12S23 



In^ 



2 /^-_Sl2 
-S23 



0(6). 



(A21) 



6. One off-shell leg 



h{D = 4 - 2e; 1, 1, 1, 1, {si2, S23, 0, 0, 0,p^}) 



Cr 



X 



{{-si2r + (-^23)-^ - - 2 Li2(i - ^; 

e V ^512 



S12S23 

2 



2Li2(l-^) 

■S23 



-Sl2\ 



"-523 



+ 0{e) . 



(A22) 



c. Two contiguous off-shell legs 



h{D = 4 - 2e; 1, 1, 1, 1, {si2, S23, 0, 0,p^,pi}) 



Cr 



■S12S23 



X 



2LMl-iiV2Ljl-iiV'°' f— 

-523/ V -523/ \—S22, 



f 0{e) . (A23) 



d. Two opposite off-shell legs 



h{D = 4 - 2e; 1, 1, 1, 1, {^12, S23, 0,^^, 0,^^}) 



Cr 



(S23S12 -vlvl) 



- (1-^12 

2Li2 ( 1 - ^"j - 2Li2 ("l - ^"j - 2Li2 (1 - 2Li2 f 1 - ^ 



■S23 



+ 



2Li2 (l-M\^x^^ (Zl21\]+0{e). 

V ■512S23/ \~S23/ -I 



(A24) 



The pole at S23S12 = p\p2 is only apparent because it is canceled by the numerator. In fact, 
setting P2P4 = ■§12523(1 — t) and performing the expansion in e we find that the integral may 
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be approximated as 



hiD = 4 - 2e; 1, 1, 1, 1, {su, ^23, 0,p^, 0,p^}) 



Cr 



S23S12 



-(2 + r) + 2 - ^ + (2 + r)(ln(-si2) + ln(-S23) - H-pD) 



+ 2Lof^)+2Lof^l+r 
V-S23/ V-S12 

where Lq,Li are defined as, 



"■523 



-S12 



J I ^ ln(r) Lo(r) + 1 



0{ey), (A25) 



(A26) 



1 — r 1 — r 

The functions Lq and Li have the property that they are finite as their denominators vanish. 



e. Three off-shell legs 







The resuh for the box integral with three off-shell legs is 

•2 2 „2 1 \ '^r 



h{D = 4 - 2e; 1, 1, 1, 1, {si2, S23, ^.pipipl)) 



[S23S12 - pipi) 



X 



+ 



^[[-si2r + [-S23) 



-P2) 



-pir- - i-pl)-') 



^((-P2)~^(-P3)^^)/(-^23)"^ + ^((-pD"^(-P4)^^)/(-^12)-^ 



- 2112(1-^1-2112^1 



2 Lio 1 



2 2 
P2P4 



■5l2'523 / 



-^12 
"■523 



(A27) 
0{e). 



f. Four off-shell legs 



In the case of four off-shell legs we have the remarkable relation 



h{D = 4 - 2e; 1, 1, 1, 1, {si2, S23,Pi,i?2)P3'P4}) = h{D = 4 - 2e; 1, 1, 1, {312523,^1^3,^2^4})) 

(A28) 

which is valid for the case when all momenta are space-like. For the continuation to other 
regions we refer to ref. 



49^. 



5. Pentagon Integrals: I5{D = 6 — 2e; 1, 1, 1, 1, 1) 



We choose as our final basis integral the finite six- dimensional pentagon, I^lD = 6 — 
2e; 1, 1, 1, 1, 1). However at one loop the results do not depend on the value of this finite 
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integral, since it always appears with a coefficient of O (e) 32l l50|. 
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